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Abstract. We study transport properties of anisotropic Heisenberg model in a 
disordered magnetic field and in the presence of dephasing due to external degrees of 
D | freedom. Without dephasing the model can display, depending on parameter values, 

?— j ■ the whole range of possible transport regimes: ideal ballistic conduction, diffusive, or 

ideal insulating behavior. We show that the presence of dephasing induces normal 
diffusive transport in a wide range of parameters. We also analyze the dependence of 
spin conductivity on the dephasing strength. In addition, by analyzing the decay of 
spin-spin correlation function we find a long-range order for finite chain sizes. All our 
results for a one-dimensional spin chain at infinite temperature can be equivalently 
rephrased for strongly- interacting disordered spinless fermions. 

> ' 1. Introduction 

cn 
. 

Theory of free fermions has been very successful in explaining properties of many 
condensed matter systems. Sometimes though, in particular in low dimensional systems, 
the picture of isolated noninteracting fermions is an oversimplifications and one has to 
ON include additional interactions. Integrable model of free fermions can be upgraded by 

taking into account several effects: (i) interaction - fermions can interact with each 



other, (ii) disorder - fermions can experience different local eigenenergies and, (iii) 
coupling - fermions can interact with external degrees of freedom. A non-perturbative 
inclusion of any of these effects greatly complicates the analysis forcing one to use various 
successful phenomenological theories [T]. There is nevertheless a desire to understand 
properties of such strongly-interacting many-body disordered systems from the first 
principles, that is starting from the governing equations of motion. 

In the present work we shall study transport properties of a one-dimensional (ID) 
antiferromagnetic Heisenberg model at infinite temperature in the presence of all the 
above mentioned effects: interaction, disorder and coupling to the environment. To our 
knowledge this is a first study of a large many-body system in such setting by using 
microscopic equations of motion. In terms of Pauli matrices the ID antiferromagnetic 
spin- 1/2 Heisenberg model reads 

ra—l n 

# = £ K^i + ^K+i + A ^ z +i) + £ ^ (i) 

3=1 3=1 
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Figure 1. (Color online) Parameter space of the Heisenberg model in the presence of 
interaction A, disorder with strength h and dephasing 7. Rigorous transport properties 
are known only on axes A = or h = 0, both with 7 = 0. In the present work we 
show that dephasing (7 7^ 0, yellow balls) induces diffusive transport. 



Magnetic field strength hj at site j is chosen according to a uniform distribution in the 
interval [— h, h]. Even though a clean Heisenberg model without disordered magnetic 
field is explicitly solvable by the Bethe ansatz and has been studied for several decades, 
its spin transport properties are far from being understood, for an overview see [2]. 
By using the Jordan- Wigner transformation it can be mapped to a system of strongly 
interacting spinless fermions, with A being the interaction strength. Therefore, the 
problem of spin transport in the Heisenberg model is equivalent to that of particle 
(charge) transport in a system of interacting fermions. In the absence of coupling 
to the environment, rigorous statements are available only for special cases without 
anisotropy (interaction), A = 0, or in the case of no disorder, h = 0, see Fig. [TJ For 
A = the model is equivalent to a system of non-interacting fermions and can be 
described in terms of single-particle states that exhibit exponential localization in the 
presence of nonzero disorder, h ^ 0, resulting in a perfect insulating behavior due to 
Anderson localization [3]. If disorder is zero the Heisenberg model displays ballistic spin 
transport in a gapless phase for A < 1 [1], while it is likely a diffusive spin conductor 
for A > 1 [HJEJITIIBIE]. This last result about diffusive transport in an integrable 
model, based mainly on numerical investigations, is somewhat unexpected and still 
lacks a deeper understanding, see though [TQl [11] for theoretical arguments supporting 
diffusive transport. As soon as interaction A and disorder h are both nonzero things 
get more complicated because one has to deal with a genuine many-body physics, 
for single-particle treatment see, e.g., (T2J [13]. There are in fact opposing claims in 
the literature on whether the localization is destroyed or preserved in the presence of 
strong interaction; some predict a phase transition from an insulating to a conducting 
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regime at high temperature [2], other numerical simulations, albeit on small systems, 
suggest localization even at an infinite temperature [T5| [16] or even ideal conducting 
behavior [T7] . 

Other non-equilibrium properties, like dynamics in various quantum quenches, has 
also been studied intensely in recent years, see [IE] and references therein. For a recent 
study in 3D see [19]. Besides its theoretical importance, Heisenberg model is also 
experimentally realized in spin-chain materials [20] . It is argued that an anomalously 
high spin conductivity of the Heisenberg model is required to explain large measured 
heat conductivity. While a detailed engineering of interactions in condensed-matter 
systems is probably out of question, experiments with cold gases in optical lattices 
could in near future lead to controlled experiments with strongly correlated many-body 
systems [21]. It is therefore important to get a better understanding of the influence 
of the interplay between interaction, disorder and external coupling on the dynamics in 
such systems. 

2. The model 

We are going to study transport by numerically solving dynamical equations of motion, 
finding a non-equilibrium stationary state (NESS) to which the system converges after 
long time. Once we calculate NESS we can evaluate various expectation values in that 
stationary non-equilibrium state. The dynamics of the system will be described by the 
Lindblad master equation [22J, 

±p = i[p,H] + £ h ^(p) + £ dc »\p). (2) 

Bath superoperator £ bath takes into account the coupling to unequal heat baths at both 
ends, inducing a non-equilibrium situation, while £ de P h represents the effect of dephasing 
caused by the coupling to some external degrees of freedom. Even though description 
in terms of Lindblad equation, originating in quantum optics, has not been employed 
very often in condensed matter physics, it has been recently used in several works 
studying spin transport in the Heisenberg model [231 I2U El [25] . While a number of 
assumptions is involved in the derivation of the Lindblad master equation from unitary 
evolution of a combined system [26], Markovian approximation can be justified because 
we are interested in the stationary state reached after long time, e.g., much longer than 
the bath correlation time. Apart from coupling to the bath the only non-unitary effect 
we take into account is dephasing. At high temperatures, which is the regime studied 
in the present work, dephasing is the dominant environmental effect. We are going 
to solve master equation using time-dependent density matrix renormalization group 
(tDMRG) method enabling us to study by an order of magnitude larger chains [TJ E] 
than is possible with other methods. Dephasing present in our simulations in fact acts 
stabilizing on the numerical method enabling us to simulate relaxation in chains of up-to 
512 spins. 
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Figure 2. (Color online) Hamiltonian part acts on all spins, Lindblad part for the 
baths acts on the boundary two spins, while the dephasing part acts independently on 
each spin. 



Let us describe various terms in the Lindblad equation (jSJ). The Lindblad bath 
superoperator £ bath affects only two border spins at the left and right end of the chain 
and involves four Lindblad operators L\ =1 2 3 4 at the left end and four Lf =1 2 3 4 at the 
right end, 

£bath(p) = ^ ( [L L,R p) L L,R t] + [L L,R pL L,R t] ) (3) 

fc 

L R, 

Lindblad operators L k ' are chosen in such way that if they would be the only term 
present in the master equation, the NESS, i.e., the state p(t — » oo), would be equal to 
the grand canonical state on two border spins. They are obtained by targeting the grand 
canonical state at an infinite temperature, p gra nd. ~ exp (— /i Lj r£ z ), S z = J2k°~h on l wo 
boundary spins, see [211 f° r more details. By choosing different potential /xl and /xr 
at the left and the right end, a NESS with a nonzero spin current is obtained. The bath 
operators are such that for /i L = /i R they induce the correct equilibrium grand canonical 
state in the bulk of an interacting chaotic spin chain [27] . Dephasing part £ de P h ; due to, 
for instance, coupling with phonons, involves only one Lindblad operator [j] L = \o~^oj 
for each spin site j, 

£ dcph (p) = ^ E {Wtcrjp, a+aj] + [a}*T, p a/,:]) , (4) 

where the sum over j runs over spin sites, a ± = a x ± i a y , and 7 is a dephasing strength. 
Dephasing (also called phase damping or phase-flip in quantum information) causes an 
exponential decay of the off-diagonal elements in the diagonal basis of a z (in fermionic 
language this corresponds to the off-diagonal decay in the number basis). It might 
be interesting to note [28J that the dynamics of operators in the Heisenberg model in 
the presence of dephasing, i.e., evolution by the Lindblad equation (j2J), can be exactly 
mapped to the evolution of pure states in a one- dimensional non-hermitian Hubbard 
model with complex on-site repulsion (iU)ntfTij±. 

I The same £ d °P h is obtained for L = ^a^. 
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Recently, the influence of dephasing on transport has been studied in short 
disordered networks modelling the light-harvesting complex [291 USB EH [32], where the 
presence of dephasing can enhance network's ability to transmit excitations. 

3. Non-equilibrium steady states 
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Figure 3. Spin and current profiles in the stationary state of a chain of n = 128 spins 
with A = 1.5, disorder strength h — 2 and dephasing 7 = 1. Bottom panel shows the 
size of magnetic field at individual sites. 

First we are going to show expectation values of some simple observables in NESS 
obtained as an infinite-time solution of the Lindblad equation (|2]). We choose bath 
potentials on left and right ends as /il,r = ±0.1 (same throughout the paper) and 
numerically solve the Lindblad equation, starting with some arbitrary initial condition. 
Details of our numerical implementation of tDMRG can be found in the Appendix. 
After sufficiently long time the state p(t) converges to a time-independent NESS. This 
"relaxation" time after which we converge to the NESS depends on the bath details, 
i.e., it is a combined property of the central system and of the coupling. It increases 
with the increasing chain size n. For each run we carefully check that the simulation 
time has been long enough so that the convergence to the NESS is indeed reached. 
We show in Fig. [3] the expectation value of local magnetization a\ and of local spin 
current jk = 2{a k a y k+l — c y k a k+1 ) in a NESS of a disordered Heisenberg model with 
n = 128 spins. One obtains a linear spin (magnetization) profile, as it is expected for a 
diffusive conductor, as well as a homogeneous spin current jk throughout the chain. Our 
numerical experience tells us that the homogeneity of the spin current can be used as a 
rather sensitive indicator of the convergence to NESS as well as of truncation errors. 
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3.1. Spin conductivity 

To determine whether the system is a normal or an anomalous conductor one has to 
study the scaling of the spin current with the system size. First, we are going to 
study clean system without any disorder, h = 0. Dephasing strength is always 7 = 1 
and three different anisotropies A = 0,0.5, and A = 1.5 are used. Keeping chemical 
potentials fixed we have calculated the spin current in the NESS for systems between 
n = 8 and n = 256 spins long. Results are shown in Fig. HJ We can see that at 
7 = 1, regardless of the value of the anisotropy, the current always scales as ~ 1/n 
with the system size, meaning that the system is diffusive (i.e., normal) conductor. The 
asymptotic value of the product j ■ n determines the coefficient of spin conductivity 
k through j = —k{{g\) — (a^))/n. For large n we have (af) — (a^) w 0.153 giving 
k = 3.79,3.20,1.53 for 7 = 1 and A = 0,0.5,1.5, respectively. With the increasing 
interaction A the conductivity decreases. For A = 1.5 and without dephasing and 
disorder, 7 = and h — 0, the coefficient of spin conductivity has been calculated in 
Refs. [7J |9] to be k = 2.30 and is therefore larger than with dephasing. 
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Figure 4. Scaling of the spin current with system's size in a clean system, h = 
and 7 = 1, at a fixed driving potential. For all anisotropies the presence of dephasing 
induces normal (diffusive) transport, signaled by the scaling j ~ 1/n (straight lines 
in the main plot). In the inset the convergence of the product j ■ n to its asymptotic 
value is shown. 



Similar results are obtained also in the presence of disorder. In Fig. [5] we show the 
scaling of the spin current with n when each site, in addition to the dephasing with 7 = 1, 
experiences also a random magnetic field hj with strength h — 2 [§. The difference in 



§ Because fluctuations between different realizations of disorder were found to be small for large n all 
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magnetizations between left and right ends is again almost independent of A and equal 
to (erf) - «) « 0.153, resulting in k = 1.63, 1.50, 1.00 for A = 0,0.5, 1.5, respectively. 
Compared to the clean case disorder decreases conductivities, however, transport is still 
diffusive. Dephasing therefore breaks any possible many-body localization present at 
7 = 0. In the inset of Fig. [5] we also show data for the case of single-particle Anderson 
localization (7 = 0, A = 0) for which the current exponentially decreases with n. 
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Figure 5. Same as Fig.2]for a disordered system, h — 2 and 7 = 1. In the inset data 
for Anderson localization in the absence of dephasing and interaction, 7 = 0, A = 0, 
is also shown. 

Dephasing therefore induces normal transport irrespective of the anisotropy and 
disorder strength. Decreasing dephasing strength to one of course has to recover known 
behavior for 7 = which is, depending on A and h (see also Fig. [I]), superconducting (for 
A < 1 and h = 0), diffusive (for A > 1 and h — 0) or insulating (for A = and h ^ 0). 
To elucidate this transition we have studied how the coefficient of spin conductivity k 
changes as the dephasing strength 7 is varied. For each value of parameters A and h we 
calculated NESS for sizes n = 64, 128 and n = 256, from which we determined n. Spin 
current in all cases scaled as ~ 1/n, the transition to asymptotics though happening at 
larger n for smaller values of 7. Dependence of k in the absence of disorder is shown 
in Fig. [61 One can see that for A = 0.5 spin conductivity increases with decreasing 
7. Divergence is scaling as k ~ I/7 which is consistent with the superconducting (i.e., 
ballistic) limit at 7 = 0. For A = 1.5 though, k converges to a finite value as 7 — > 0, 
in accordance with the normal transport in the Heisenberg model at A > 1 without 
dephasing. 

data shown are for a single realization of disorder. 
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Figure 6. Dependence of the spin conductivity k on the dephasing strength 7 in the 
absence of magnetic field. 

In the presence of disorder the limit 7—7-0 corresponds to Anderson localization if 
A = while for nonzero A things are not so clear. From our results seen in Fig. [7] we 
see that for small 7 spin conductivity decreases, however, the limit value ^(7 = 0) is 
rather difficult to predict. Spin conductivity reaches a maximum at some intermediate 
value of 7 and again decreases for large 7. Note that there is no qualitative difference 
between A = 0.5 and A = 1.5. A nontrivial maximum of conductivity is found also 
when studying dephasing in photosynthetic complex j29l [30j EH |32] . 
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Figure 7. Dependence of the spin conductivity k on the dephasing strength 7 in the 
presence of disordered magnetic field with strength h = 2. 
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Figure 8. (Color online) Absolute value of the spin-spin correlation function C(i,j) 
for NESS of a chain with n = 128 spins and h = 0, 7 = 1. Color code is for 
log 10 \C(i, j)\ (left panel). Right panel: cross-section along the diagonal, C(r) = 
C(n/2 - r/2,n/2 + r/2); data are for n = 64,128,256 and n = 512 (left to right). 



3.2. Correlation function 

We have also calculated the spin-spin correlation function for non-equilibrium steady 
states from Figs. H] and [5j Results are similar for all values of A and we show only the 
case with A = 0.0, for which tDMRG is the most efficient. Spin correlation function 
C(i,j) = (cr^o-j) — (cTj Z )(cr|), shown in Figs]8] and [9j in all cases has a plateau. This 
indicates the presence of a long-range order. Scaling analysis shows that the plateau 
in C(i,j) scales with the system size as ~ 1/n and with the potential difference as 
~ (A/i) 2 = (//l — /Ur) 2 - The same scaling is obtained also for A 7^ (data not shown), 
however, because of less favorable scaling of the entanglement with n we could reliably 
calculate C(i,j) only for sizes n ~ 128. Note that decreasing values of the correlation 
function for increasing size make a precise calculation of C(i,j) rather demanding. 
Scaling of the correlation function plateau C ~ (/il — /^r) 2 / n ~ j ■ (//l — A*r) means that 
the long-range order is a purely non-equilibrium phenomenon. In equilibrium, when the 
spin current is zero, j = 0, it goes to zero. In addition, correlations go to zero also 
in the thermodynamic limit n — > 00, as one would expect for a true normal (diffusive) 
conductor. Similar scaling of the correlation plateau has been recently observed in a 
non-equilibrium phase transition in the XY spin chain in a homogeneous magnetic field, 
studied by exactly solving the corresponding master equation [33| l3l]. 

4. Conclusion 

We have studied spin transport in the spin-1/2 Heisenberg model at infinite temperature 
in the presence of disorder and dephasing. By solving master equation for non- 
equilibrium driving potential we numerically obtained a non-equilibrium steady state 
and calculated expectation values of different observables. Studying the scaling of the 
spin current with system's size for systems of up-to 256 spins we found that the dephasing 
always induces diffusive transport regardless of the anisotropy or disorder strength. The 
size of spin conductivity coefficient has a non-trivial maximum at intermediate values 
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Figure 9. (Color online) The same as in Fig. El but for a disordered model with h = 2. 



of the dephasing strength for a disordered model while it monotonically decreases with 
dephasing in the absence of disorder. Studying the spin-spin correlation function in a 
non-equilibrium steady state we have found the evidence for a long-range order at all 
anisotropies. The value of the correlation plateau scales as ~ — /iR,) 2 /n, meaning 
that the correlations disappear in the thermodynamic limit as well as in the equilibrium 
situation. It would be interesting to study more in detail how the transition from 
diffusive transport to either superconducting, diffusive or insulating behavior happens as 
the dephasing strength is decreased to zero. More studies would be also needed to clarify 
the situation with a non-zero anisotropy and without dephasing, which corresponds to 
strongly interacting fermions. 
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Appendix 

Here we briefly describe our implementation of tDMRG, for more details see |36j . 
Appendix A.l. Matrix product formulation 

Any state p from a Hilbert space H which is a tensor product of n local spaces of 
dimension d, W = "H^™., dim(7/i oc .) = d, can be written in a matrix product operator 
(MPO) form, where expansion coefficients are expressed in terms of product of matrices. 
If we denote basis states of T-L\ oc . at j-th site by crj, v — 0, . . . , d — 1, we can write 

\p) = Y,JM?M? ■ ~M?x\o?o? ■ ••<"). (A.l) 

Ui 

For each site j we have d matrices Mj J , each of dimension D x D while x is some 
arbitrary D dimensional vector [jj, in our case its components are X{ = 5^1. When 

|| Frequently one uses a trace instead of projection on x. For large D two formulations are the same. 
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solving the Lindblad equation (j2J), \p) represents a density operator which is viewed 
as a member of a Hilbert space of operators, with local basis being Pauli matrices 
cr° = (j x , a 1 = a 7 , a 2 = o~ z ,a 3 = 1. For a given \p) the choice of matrices Mj J is not 
unique, however, as we will see, there is a preferred choice corresponding to Schmidt 
decomposition of \p). 

Appendix A. 2. Evolution 

The above setting is used to solve master equation p = Cp, where £ is some linear 
operator. In our simulation £ is just the operator corresponding to the right side of 
the Lindblad equation (j2J). Solution can be formally written as p(t) = ]l ex P (CAt) p(0) 
in terms of propagators exp (£At) for a small time step At; in our simulations we use 
At = 0.05. Each small-step propagator is then decomposed according to the Trotter- 
Suzuki formula. Writing £ = A + B as a sum of two parts, each of which is a sum of 
mutually commuting terms, we have 

exp (£At + 0(At k+1 )) « exp ( ai AAt) exp (hBAt) exp (a 2 AAt) ■ ■ ■ exp (b k BAt). (A.2) 

We use a fourth order {k — 4) decomposition [35] with a\ — s/2, b\ — s, a 2 — 
(l-s)/2, b 2 = (l-2s), a 3 = (l-s)/2, b 3 = s, a 4 = s/2, b 4 = 0, where s = 1/(2-^2). 
In our Lindblad equation we have only nearest neighbor unitary terms and at most 
two-qubit nearest neighbor dissipative terms. We put into A all even bonds, into B all 
odd bonds, while single qubits terms are equally split between the two terms. A and 
B are therefore sums of commuting two-qubit terms. Basic operation we have to do 
is then exp (Ajj+ie), where e oc At. For time-independent Ajj+i this is a fixed linear 
operator on the Hilbert space of operators. Its matrix representation A/j+i in the basis 
joy 7 ) can be calculated in advance, 

[Aj, j+1 ]^v j+1 ,v jVi+1 = tr (o-ytf+i exp (Ajj+^a^aj^ 1 )/^ (A.3) 

Ajj +1 is orthogonal for unitary exp (A, J+1 e). 

Every bipartite state can be written in the form of the Schmidt decomposition. 
Splitting a system into the first j sites (1, . . . , j) and the rest, we can write 

|p) = E^' ) K ' ) )®K°' ) ), (a.4) 

a 

with an explicit form of orthogonal states \w^^) and \w^^}, 

l«4 w > = E xpW? ■ ■ ■ ^W-^K ■ • • of) 

f3,Ul,...,Vj Aq 

= E MS 1 ■ ■ ■ M?]aj> *P\°&L ■ ■ ■ <">■ (A.5) 

In terms of matrices Mp the orthogonality of \w^^) is reflected in the condition 

E {M; j } t diag({A (j - 1) } 2 )Mf = diag({A (j) } 2 ), (A.6) 
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for all j, while for |u>„ ) 

^MpiMpY = 1, (A.7) 

must hold. Coefficients A^ are called the Schmidt coefficients. How well the tDMRG 
method works crucially depends on how fast the ordered A^ decrease with a. In general 
time evolution can not be done in an exact way unless we increase the dimension of 
matrices M- j exponentially with time. To keep dimensions fixed we have to make a 
truncation. This is done in an optimal way (in terms of probabilities (A^) 2 ) if we make 
sure that we drop states corresponding to the smallest Schmidt coefficients. To do this 
we first calculate the action of on the matrix product operator (MPO) form (1A.1I) ; 

it transforms a product of two matrices Mp and iWj^ 1 at sites j and j + 1 into a single 
bigger matrix O, 

e» ja ,» J+1 p = xt 1] E [Ajj+xU^,^ [m?m;^u. (a.s) 

The reason to multiply it with A^ -1 -* is to recover Schmidt decomposition after doing a 
singular value decomposition on 0, decomposing it as G = UdV, with unitary U and V 
and diagonal d. It turns out that if we prescribe new matrices M and new coefficients 
A^ as 



7 




J J a, 7 


A(i) 

= U 7 

u »ia,i (i _i) 
Aa 


J J 7,a 





(A.9) 



and if exp (Aj j+ie) is unitary, orthogonality conditions (IA.6|IA.7[) are preserved. The 
whole Trotter- Suzuki step ( IA.2j) . advancing in time by At, takes 0(kn(dD) 3 ) operations. 

Appendix A.S. Reorthogonalization of MPO 

One Trotter- Suzuki step is optimal if all transformations preserve Schmidt 
decomposition, i.e., if they are unitary. Coupling to the bath and dephasing obviously 
violate this condition. Because the bath part £ bath destroys unitarity only at the two 
boundaries we included it in A and B terms. Dephasing part though, acting on each 
spin, would affect unitarity also in the bulk of the chain. As a matter of numerical 
convenience we have not included it in the Trotter-Suzuki decomposed part flA.2|) but 
have instead applied it only at the end of each Trotter-Suzuki step as exp (£ deph At). 
This corresponds to dephasing acting in a kicked manner and not in a continuous way 
and has no physical consequences for our results. 

After each step of length At, we are left with a MPO \p) whose matrices do 
not correspond to Schmidt decomposition anymore (due to non-unitary propagators in 
between). Before proceeding to the next step we reorthogonalize matrices Ma 3 in order 
to recover the optimality of the Schmidt decomposition. This takes of order 0(ndD 3 ) 
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steps and is done in the following way [3T| I3"5] . First, we recursively construct matrices 
V]/ and Vf as 

yf = 1 Z(M?yV?L 1 M?, j = l,...,n-l (A.10) 
V£ 1 = J2M?V*(M?)\ j=n,...,2, (A.11) 

starting with the matrix [V^]^/ = x^xf and [V*] ki i — Xfcx\ at the boundaries. If 
conditions ( 1A.6IIA.7I) are satisfied the matrices Vj'Vj 1 are diagonal. In general this 
is not the case so we have to rotate each Mp in order to recover orthogonality. To 
do this we calculate a square root of a non-negative Vj" and diagonalize the matrix 
^ = \A^^ R \A7" = U dW, in terms of diagonal d and unitary U. Finally, we calculate 
unitary d = d^W (Vf) 1 ' 2 and its inverse Gj 1 = (V^'^Ud 1 / 2 for j = 1, . . . , n - 1. 
We also set G = G n = t. New matrices Mp in MPO (lA.lj) . respecting Schmidt 
decomposition, are then obtained by rotations 

Mp = G^MpGj 1 , j = 1, ... ,n. (A.12) 
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